#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <pthread.h>

// Número de iterações = 10^5
#define M 100000

// Número de threads
#define N_THREADS 8

pthread_t threads[N_THREADS];
pthread_mutex_t mutex_sum;

double S, E, r, sigma, T, trials[M], sum = 0;
unsigned short xsubi[3];

void *blackScholes(void *offset){

	double t, max;

	// Divisão do número de vezes (intervalo) que cada thread vai rodar
	// Primeira thread: for 0 to 12499; segunda thread: for 12500 to 24999; e assim por diante
	int index = *((int *) offset);
	int div = M / N_THREADS;
	int from = div * index;
	int to = from + div;
	int i;
	int partial_sum = 0;

	// Cálculo do valores do vetor "trials" segundo a equação de Black Scholes
	for(i = from; i < to; i++){
																	// erand48 = geração de números aleatórios [0.0, 1.0)
		t = S * exp((r - (sigma * sigma) / 2.0) * T + sigma * sqrt(T) * erand48(xsubi));

		max = 0.0;
		if(t - E > 0.0)
			max = t - E;

		trials[i] = exp(-r * T) * max;

		// Para que não seja necessário utilizar o lock do mutex "mutex_sum" em todas as iterações do loop para cada thread, é feita uma soma parcial
		// Dessa forma, como as variáveis locais não são compartilhadas entre threads, só é necessário um lock para atualizar a soma final
		partial_sum += trials[i];
	}

	// Para que somente uma thread accesse a variável global "sum", é utilizado o mutex "mutex_sum". Assim, a variável "sum" é atualizada corretamente e de forma eficiente
	pthread_mutex_lock(&mutex_sum);
	sum += partial_sum;
	pthread_mutex_unlock(&mutex_sum);

	pthread_exit(NULL);
}

int main(int argc, char *argv[]){

	double mean, stddev, confwidth, confmin, confmax;
	int i, rc;

	scanf(" %lf %lf %lf %lf %lf", &S, &E, &r, &sigma, &T);

	pthread_mutex_init(&mutex_sum, NULL);

	for(i = 0; i < N_THREADS; i++){
		// Criação das threads
		rc = pthread_create(&threads[i], NULL, blackScholes, &i);

		if(rc){
			printf("ERROR: %d\n", rc);
			return -1;
		}
	}

	// Junção das threads quando completarem sua execução (em ordem de criação)
	for(i = 0; i < N_THREADS; i++)
		pthread_join(threads[i], NULL);
	
	// Cálculo da média
	mean = sum / (double) M;
	sum = 0.0;

	// Somatório para futuro cálculo do desvio padrão
	for(i = 0; i < M; i++)
		sum += (trials[i] - mean) * (trials[i] - mean);

	// Desvio padrão
	stddev = sqrt(sum / ((double) M));

	// Intervalo de confiança
	confwidth = 1.96 * stddev / ((double) sqrt(M));
	confmin = mean - confwidth;
	confmax = mean + confwidth;

	printf("S\t\t%.0lf\nE\t\t%.0lf\nr\t\t%.0lf\nsigma\t%.0lf\nT\t\t%.0lf\nM\t\t%d\nConfidence interval: (%.6lf, %.6lf)", S, E, r, sigma, T, M, confmin, confmax);

	pthread_mutex_destroy(&mutex_sum);
	return 0;
}
